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1.0 Introduction 

This report concerns the use of multiquadric functions to 
approximate scattered data. Here we deal with functions of two 
independent variables, but the methods are easily extendible to 
arbitrary dimensions, and we expect that many of the conclusions 
will carry over. 

The impetus behind our investigation is that of obtaining 
surface approximations that are efficient in subsequent 
applications. That is, we consider it to be acceptable to expend 
considerable computational resources to obtain the approximation 
in a preprocessing step. Once obtained, the approximation should 
be able to be evaluated fairly efficiently such as when it is to 
be used numerous times in an application program. 

The scattered data approximation problem is easily described 
and occurs frequently in many branches of science. The problem 
occurs in any discipline where measurements are taken at 
irregularly spaced values of two or more independent variables, 
and is especially prevalent in environmental sciences. We will 
suppose that triples of data, (Xj,yj,Zj), j=l, ...f N are given, 
assumed to be measurements (perhaps with error) of an underlying 
function z=f(x,y). The function f is to be approximated by a 
function F(x,y) from the given data. A recent survey of such 
methods is given in [FN91] . 

Multiquadric functions were introduced for interpolation of 
scattered data by Hardy [HA71]; also see [HA92] for a historical 
survey and many references. The method is one of a class of 
methods known now as "radial basis function methods" that 
includes other attractive schemes such as thin plate splines 
[HD71, DU76, and others]. The basic idea of such methods is 
quite simple, and we describe it in some generality; for 
purposes of being definite it is pertinent to note that for the 
multiquadric method the radial function is h(d) = V(d^+r^) . In 
general, suppose a function of one variable, h(d), where d 
denotes distance, is given. 

For interpolation (that is, exact matching of the given 
data), a basis function, Bj(x,y) = h(dj) is associated with each 
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data point. Here dj = V ( (x-Xj ) (y-yj ) ^) , distance from (x,y) to 
(Xj,yj). Thus each basis function is a translate of the radial 
function, h. The approximation is a linear combination of the 
basis functions, along with some polynomial terms that may be 
necessary in some cases, or may be used to assure that the 
approximation method has polynomial precision. Thus, 

N M 

(1) F(x,y) = S ajBj(x,y) + S bjqj(x,y) 

j=l j=l 

where {qj} is a set of M polynomials forming a basis for 
polynomials of degree <m. The coefficients aj and bj are 
determined by the linear system of equations prescribing 
interpolation of the data, and exactness for polynomials of 
degree <m: 

N M 

S ajBj(Xi,yi) + S bjqj(Xi,yi) = 2^, i=l, ..., N 
j=l j=l 

( 2 ) 

N 

2 ajqi(Xj,yj) = 0, i=l, ..., M . 

j=l 

For multiquadric basis functions, this system of equations is 
known to have a unique solution for distinct (Xj,yj) data (see, 
for example, [MI86]); while m may be taken as zero (no 
polynomial terms) , the theory indicates that a constant term 
should be included, and we have done so in all our work. If 
higher degree polynomial precision is desired, inclusion of those 
terms imposes no particular burden. 

While interpolation theory is important and indicates 
something about the suitability of the class of functions for 
approximation purposes, our emphasis here is on least squares 
approximation. This implies using fewer basis functions than 
there are data points. In analogy with univariate cubic splines, 
it is convenient to refer to the points are which the radial 
basis functions are centered as "knots", as was done in [MF92], 
and we do so here. If a set of knot points, (U)^/Vj.) , k=l,...,K , 
with K<N have been specified, then the problem of fitting a 
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multiquadric function by least squares is similar to that of 
solving the system of equations corresponding to those above in 
the least squares sense. We give the details. Now, let Bj^(x,y) 
denote the radial basis function associated with the point 
(uj.,Vk) / Bk(x,y)=V( (x-Uk) ^+(y-Vk) ^) . The system of equations, 
specialized for our case, is now of the form 

K M 

2 ( Xi , y i ) + 2 b^q^ (Xj^,yj^) — i=l, . . . , N 

k=l j=l 



S ak — 0 . 
k=l 

There is a question of how to treat the last equation, which 
guarantees polynomial precision. In [FC92] the corresponding 
constraint equations were imposed exactly, rather than 
approximately because of physical considerations. While there is 
not the corresponding physical situation here, we have also 
imposed the last equation as a constraint. This constraint can 
be used to reduce the size of the system by solving for aj^ in 
terms of the other ak and substituting into the first set of 
equations . 

If the knot points are a subset of the data points, then the 
same theory that assures a unique solution of the system for 
interpolation also guarantees a unique solution of the least 
squares problem. When the knot locations may differ from the 
data points, the problem of whether the coefficient matrix is of 
full rank or not is unknown to us, although we feel certain that 
the matrix is of full rank when the knot points are distinct, and 
have encountered no situations that indicate otherwise. 

In our implementation of the algorithms described in 
subsequent sections we have used a QRP' decomposition of the 
coefficient matrix to solve the least squares problem. This 
provides a stable and efficient means for solution of the problem 
with an indication if a matrix of less than full rank is 
encountered. 

In order to test the algorithms we have used a number of 
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data sets. Several of these are based on previously published 
and widely available (x,y) data sets and parent functions. We 
have also used a few less readily available data sets that we are 
willing to share with anyone interested in obtaining them. Table 
1 gives a summary of most of the data set. 



n.m This refers to point set n and function m from [FR82], 

for n=l, 2, and 3, and m=l, ..., 6. n=l is 25 points, 
n=2 is 33 points, n=3 is 100 points. n=4 refers to the 
200 point data set used in [MF92]. m=l is the humps 
and dip function, m=2 is the cliff, m=3 is the saddle, 
m=4 is the gentle hill, m=5 is the steep hill, and m=6 
is the sphere. In addition, m=7 refers to the curved 
valley function from [NI78]. 

GT This refers to the thinned glacier data consisting of 

678 points, with certain contour lines removed, from 
[MF92] . 

GL This refers to the thinned glacier data consisting of 

873 points. 

HF This is the data set from [MF92] generated to be 

approximately proportional to curvature, consisting of 
500 points. 



Table 1: Data sets used extensively in tests 

Section 2 deals with a "greedy" algorithm for determining 
the location of a reasonable set of knots for approximation of 
given data by a least squares multiquadric function. Some 
experiences with the method are given. In Section 3 we expand 
the algorithm to include the knot locations and the parameter 
value of the method as part of the optimization process. Some 
results and observations about the process are made, with the 
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optimized value of the parameter r being of interest. The 
occurrence of near multiple knots is a particularly interesting 
phenomenon. In section 4 we further extend the algorithm to 
include variable parameter values at the knots. The optimized rj^ 
values and the near multiple knots are again of special interest. 
Finally, in Section 5 we discuss some ideas for further 
exploration of least squares multiquadric approximations. 

2.0 An Adaptive Method for Knot Selection 

This section describes a greedy method for the selection of 
knot locations for fitting surfaces to scattered data using a 
least squares multiquadric function. As noted in the 
introduction, the use of fixed knots and parameter value with the 
multiquadric function results in a linear system to be solved in 
the least squares sense. These are solved using the QRP ' 
decomposition. The algorithm was implemented in Matlab^, giving 
access to powerful matrix-vector notation that simplifies many 
aspects of the implementation. In addition, Matlab allows easy 
interactive intervention in the process, with tabular and 
graphical information being made available as the computations 
proceeds. While an efficient implementation would also provide 
for updating the QRP' decomposition as more knots are added, we 
have not done this in the experimental program since our 
computational resources were sufficient to make it unnecessary. 

2 . 1 The Algorithm 

The algorithm proceeds as follows, with the necessary input 
being obtained by interrogation of the user. The description 
given starts after all input has been obtained. 

a) The initial step is to obtain the least squares fit by a 
constant function, the average of the data values. The two 
data points having maximum positive and maximum negative 
error are taken to be the first two knots, (u^,V 2 ) and 
(U 2 ,V 2 ). The knot counter K is set to 2. 

1 MathWorks, 24 Prime Park Way, Natick. MA 01760 
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b) The least squares multiquadric fit with K knots is obtained. 
The residuals are computed along with their rms value, the 
approximation is evaluated on a grid of points, a smoothness 
measure approximately equal to the value of the thin plate 
functional over the region is computed, and if the underlying 
function is known the rms error on the grid is computed. These 
values are then output and a perspective plot of the 
approximating surface is given. 

c) The maximum absolute value of the residuals is found and the 
location of this residual, subject to the minimum knot 
separation value, is taken to be the next knot location 
(Uj^,Vj^). At this point the algorithm proceeds to step b 
unless the maximum number of knot locations to be computed 
has been reached. 

At the termination of the program, the user can restart the 
process with any of the parameters changed, with any number of 
knot point locations, up to the total number that have been 
computed. Hard copy plots of the surfaces and tabular output can 
be obtained. 

2.2 Some Results 

One of the interesting aspects of the multiquadric method 
concerns appropriate choice of the parameter, r. Initial advice 
was to specify the value in terms of approximate data point 
separation [HA71,FR82], although even in [FR82] it was clear that 
the best value was dependent on the ordinate data as well. More 
recent work [TA85, CF91] has shown this to be the case and an 
algorithm for a ”good'' value was given in [CF91] . 

While no algorithm was implemented to obtain the best r 
value for fixed knots found by the adaptive method above, the 
flexibility of the implementation allowed for some interactive 
experimentation along these lines. In most cases investigated, 
it was found that the value of r used in the process of selecting 
the knot locations also was very close to the "best” value (that 
is the one that minimized the rms error of the residuals) for 
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that particular set of knot locations. Exceptions were when the 
number of knots was quite small (5w8), in which case the 
multiquadric method shows a striking affinity for best fits with 
r very close to zero. Of course, most surfaces with any 
complexity cannot be fit well with so few knots. Apparently the 
adaptive knot selection process is quit dependent on the r value 
used, at least enough to rule out significant improvement by 
changing r once the knots have been selected. 

While a reasonable a priori choice of the parameter r in 
this context can be made, the value of the best r is still an 
open question, and is not likely to be resolved anytime soon. As 
is pointed out by [CF91], the parameter can be used somewhat like 
a tension parameter (small values correspond to "tighter" 
surfaces) , and consequently surfaces that involve steep gradients 
will be approximated with less overshoot by selecting a small 
value of r. The tension effects are limited compared with the 
results that can be obtained using thin plate splines with 
tension (see [FR85]). Other factors enter into the selection of 
the r value, however, since small values also lead to rapid 
changes in the gradient which may be undesirable. 

One of the parameters in the knot selection process is the 
minimum separation between knots. It has been found that there 
is often an improvement by requiring some moderate separation 
between knots, for example imposing a minimum separation of .1 
or .2 for 20 knots on the [0,1]^ for point set 3. This tends to 
distribute knots more uniformly throughout the region, even when 
there are clumps of data. For comparison purposes, the rms 
errors (rmse) at the data points and over a 20x20 grid were 
computed and are given in Table 2 for several data sets. All 3.m 
examples were with r=0.3 and 20 knot points, while the HF data 
set used r=0.2 and 50 knot points. 
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data set 


minsep 


rmse (data) 


rmse (grid) 


3.1 


0.0 


0.0215 


0.0231 




0.1 


0.0262 


0.0281 




0.2 


0.0210 


0.0241 


3.2 


0.0 


0.0122 


0.0132 




0.1 


0.0162 


0.0168 




0.2 


0.0116 


0.0133 


3.3 


0.0 


0.0020 


0.0023 




0.1 


0.0014 


0.0018 




0.2 


0.0016 


0.0019 


3.6 


0.0 


0.0033 


0.0038 




0.1 


0.0026 


0.0031 




0.2 


0.0053 


0.0051 


HF 


0.0 


0.0031 


0.0040 




0.05 


0.0030 


0.0037 




0.1 


0.0032 


0.0045 



Table 2: rms errors for various separation distances 

In the case given in [MF92] where the data was specifically 
generated to reflect the curvature of the underlying surface, the 
knots computed by this algorithm tend to be gathered in regions 
where the density of data points is greatest. Figure 1 gives the 
results in one case. It shows the data points and the subset 
selected as knot points by the greedy algorithm, along with the 
contours of the parent surface, in part a. Here the minimum 
separation distance of 0.05 was imposed, resulting in a more 
regular distribution than when a zero separation distance is 
imposed. In part b the surface from which the data was sampled 
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is shown. This function is used in later examples (function 1 
from the table) ; the viewpoint is from the right center field. 
In part c the surface shown is that constructed by least squares 
fit using the knot points in part a. Part d shows the contours 
of the approximating surface. Part a can be directly compared 
with Figure 3 in [MF92], and it is seen that the distribution is 
different, and in particular does not have the nice spacing of 
that in [MF91] . Qualitatively the knot locations given here do 
reflect the density of the data, however. 

The greedy algorithm given here appears to be potentially 
useful for many problems where data subject to error is available 
and the surface must be apprc <imated using an approximation that 
is computationally as efficient as possible. A problem which we 
have considered, but which needs additional attention, is that of 
when enough knots have been generated so that the behavior of the 
underlying surface is captured without undue influence by the 
errors in the data. For now, this is mostly an unexplored idea, 
and we have more to say about it in Section 5. 

3.0 Variable Knot Locations and Multiguadric Parameter 

While the adaptive method discussed in the previous section 
seemed to perform reasonably well, it was felt prudent to check 
the performance of the scheme against one which considered the 
knot locations, along with the parameter value, r, to be 
variables over which the minimization of rms errors at the data 
points was achieved. The function to be minimized in this case 
is the same as before, but here we will state it explicitly 
rather than in the implied form where the least squares solution 
was that of the overdetermined system (3) . The minimization 
problem is 

N K 

(4) min 2 [Zj^ - 2 ^kB 3 ^(Xj^,yj^) - c]^ 
i=l k=l 

where the minimization is taken over all (Uj^/V]^), r, the aj^, and 
c (with the last equation of (3) imposed as a constraint) . As a 
practical matter, for each given knot configuration and r value. 
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the least squares solution of (3) computed as a step toward (4) . 
This results in the solution of a simpler, but equivalent problem 
since 2K parameters are eliminated from (4) by imposing the 
condition that the values of the aj and c be always taken as 
obtained from the least squares solution of (3). Hence, our 
final process is more properly written as 

N K 

(5) min min S - S - c]^ 
i=l k=l 

where the inner minimization is over the aj and c (least squares 
solution of (3), and the outer minimization is over the knot 
locations and the value of the parameter r. The global minimum 
of each of the two problems are clearly the same. Eg. (5) is the 
more restrictive, but any minimum of (4) is a local minimum of 
(5) , else a better solution is attainable for (4) . This does not 
imply that the iterative methods employed to solve (5) would work 
equally well, nor find the same local minima, when applied to 
(4) . 

When knot locations are allowed to differ from data 
locations, the guarantee of full rank of the coefficient matrix 
conferred on the system by interpolation theory no longer holds. 
As we noted in the Introduction, this has not posed any problems 
in our computations. 

3.1 Optimization Algorithms and Initial Guesses 

We have used two different nonlinear optimization schemes, 
both implemented and available as part of Matlab. One is the 
procedure FMINS that is based on a simplex procedure [W085]. The 
other is LEASTSQ that is based on the Levinberg-Marquardt 
procedure. Both of the routines appear to reliably find good 
local minima that are qualitatively similar, although LEASTSQ 
often finds a somewhat smaller rms residual and we have used it 
for most of the results given here. 

The initial guess has a strong influence on the solution 
obtained by any nonlinear optimization program. Except for a few 
experiments, we have used the results of the greedy algorithm in 
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the previous section, with a somewhat judicious guess at the 
value of r, as the initial values for the nonlinear optimization. 

3.2 Some Results 

One of the values of interest is the optimized value of the 
parameter r. For function 1 the usual values tended to be around 
0.1 to 0.2, although in some cases values outside that range were 
obtained; the smallest rms errors were obtained in that range. 
For function 2 much smaller values were obtained, generally in 
the range less than 0.05. For function 3, values in the range 
0.20 to 0.30 were prevalent. For function 6 the value obtained 
in the one computation we carried out was more than 10. It is 
tempting to try to compare our results with the best values found 
for interpolation by [CF91], and with their formula for 
approximating the best value. For the moment we can say that for 
the most part the data do not seem contradictory, although for 
function 3 our values are somewhat smaller. For function 6, the 
value we obtained was in line with computational experience in 
[CF91] in that the value is quite large. 

One very interesting aspect of the results of computing 
local minima of (5) is that, with the exception (and then not 
always) of computations involving fewer than 10 knots, the 
results involved near repeated knots, sometimes several different 
pairs with 20 or more knots, and sometimes triples of closely 
spaced knots. Because of the nonzero convergence tolerance for 
the optimization routine, by "near repeated” knots, we mean those 
that are within a distance consistent with the convergence 
tolerance. In some cases there were also other knots within 
distances of 0.02 or 0.03 for data in the unit square. 

The occurrence of near multiple knots suggests that the 
method is trying to adapt to some behavior of the surface that 
cannot be approximated locally by a single multiquadric basis 
function. The behavior of a linear combination of multiquadric 
functions at points far away is essentially the same as a single 
multiquadric. Because of the local extremum of the multiquadric 
function near the knot point, it was not immediately clear what 
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could be achieved by a linear combination of multiquadrics at 
nearby knot points. Because of this, an investigation of the 
behavior of the surface defined by terms in the approximation (1) 
corresponding to the near repeated knots was undertaken, and in 
particular, comparison with the surface defined by a single 
multiquadric having the knot point at the average of the repeated 
knots, with coefficient equal to the sum of the coefficients for 
the repeated knots. Far away, the behavior of the composite 
function must be, and is seen to be, essentially similar to a 
single multiquadric. In the vicinity of the knots (and not 
necessarily just between them) the behavior can be very 
different. 

Near multiple knots result in the coefficient matrix being 
poorly conditioned, which also allows for the possibility of 
large coefficients in the least squares solution of (3) . We are 
unable to deduce for certain whether the closeness of knots is 
required in order for the coefficients to become large, or 
whether the closeness is required to obtain the required behavior 
in other ways. In one case we looked at, the knots are within 
0.0035 of each other, the magnitude of the coefficients is on the 
order of 1125, and the condition number of the matrix is larger 
than 10^, some four orders of magnitude larger than needed for 
the magnitude of the coefficients since the data is on the order 
of one. 

It seems to be true that the most deviant behavior of the 
sum of the near multiple terms occurs when the sum of the 
coefficients for the nearby knots is relatively close to zero. 
As an example showing quite different behavior of the sum of the 
terms for two nearby knots from that of the average term, see 
Figure 2. Parts a and b show perspective plots of the two 
surfaces, while parts c and d show contours of the same two 
surfaces. The deviation is striking and make it seem reasonable 
that in order to capture local behavior, multiple knots are 
necessary since local behavior cannot be affected by basis 
functions that are associate with far away knots, and each basis 
function itself is locally a hyperboloid (of one sheet - no 
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saddles) in shape. 

Finally, we give the results of optimizing on the knot 
locations and the value of r for the 50 knot approximation 
corresponding to Figure 1. The results are shown Figure 3, and 
reveal that while many knots have moved from their initial 
positions, there density still reflects the same general pattern. 
Noteworthy is the fact that there is only one cluster of near 
repeated knots, those being the three at about (.44,. 78), near 
the dip in the surface. Those three are clustered within a 
distance of less than 0.01, while there is another knot within a 
distance of less than 0.025. Finally we note that the rms error 
for the surface was improved from 0.0037 to 0.00023 by the 
optimization process; r changed from the initial 0.2 to 0.2389. 

4.0 Variable Knot Locations Each with Variable Parameter 

The computational experience gained in the variable knot 
case, and especially the near repeated knot phenomenon lead us to 
consider whether or not the near multiple knots were occurring 
because a single value of the parameter at all knot locations was 
not necessarily appropriate. Thus, we modified the algorithm to 
allow for an independent r value, rj^, to be associated with each 
knot. The implications for the rank of the system is again not 
known. It is, however, easy to find examples of different 
parameter values that lead to singular systems in the 
interpolation problem. We believe the least squares problem is 
more robust, and have found no troublesome cases during our 
investigation . 

4.1 Some Results 

We soon discovered that the use of variable parameter values 
did not alleviate the problem of near multiple knots in the 
optimized approximation. It is interesting that when multiple 
knots occur, the parameter values for those knots are invariably 
very close to having the same value. In fact, our limited 
experience seems to indicate that most knots tend to have similar 
values of the parameter, although there are generally a few that 
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take on smaller values than elsewhere. 

Once again, the behavior of the surface in the vicinity of 
near multiple knots often reflects behavior that cannot be taken 
on by a single multiquadric. As an example of a different kind 
of behavior than illustrated by Figure 2, note that Figure 4 
shows another case where the near double knot results in a 
surface that resembles a quadric with a dimple in it. The rj^ 
values for the two knots are essentially the same. 

For comparison purposes between the greedy algorithm, the 
variable knot and parameter algorithm, and the variable knot each 
with variable parameter value, we look at a case with a few 
knots. In Figure 5 we give the results of the greedy algorithm 
for the data set 3-1 with 12 knots and an initial knot separation 
of 0.2. Parts a-d are, respectively, the point and knot set, the 
parent surface, the approximating surface, and contours of the 
approximating surface. In Figure 6 the results of the variable 
knot and parameter algorithm are given; the initial values were 
those resulting in Figure 5. Parts a-d are, respectively, the 
point set, the parent surface, the approximating surface, and the 
initial and final knots locations. The improvement is clear. In 
Figure 7 we see the results of the variable knot, each with 
variable parameter value. The parts of the figure are the same 
as for Figure 6. Here the improvement is even more spectacular. 
The values of the multiquadric parameter and the rms errors at 
the data points and on the grid are given in Table 3 . 



Algorithm 



r value (s) rmse(data) 



rmse (grid) 



Greedy 

Var kts, var r 
Var kts, var rj^ 



0.3 

0.158 

0 . 12 - 0.66 



0.0320 

0.0101 

0.0012 



0.0341 

0.0119 

0.0023 



Table 3: Data set 3-1 with 12 knots, initial separation of 0.2 

The improvement in the rms errors with variable knots is 
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significant for this particular data set. In situations where 
the number of knots is sufficient to give a reasonable 
approximation we find the typical improvement in rms errors is 
about by a factor of 3-10 when variable knots and parameter 
values is allowed, and another factor of 3-10 when the parameter 
values are allowed to vary. This is, however, highly dependent 
on the parent surface, for example, the cliff surface 
approximations are improved by smaller factors, while the saddle 
surface approximations tend to be at the upper end of the scale. 

It is interesting to compare the results on this particular 
example with those of [CF91] with interpolation to the same data. 
There it was found that the •' 'est'' value of the parameter r (that 
being about 0.33) lead an approximation (which is the sum of 100 
multiquadric terms) which has an rms error of 0.0026. It is 
startling to see that the 12 term approximation derived using 
variable knots each with variable r has smaller rms error. We 
have not followed this line of investigation very far, but 
function 2 (cliff) is also approximated well using relatively few 
knots . 

It appears that the use of variable knots can give a greatly 
improved approximation when using multiguadric functions with a 
fixed number of knots. When variable parameter values are 
allowed the complexity of evaluating the approximation is 
essentially unchanged and seems to be a worthwhile improvement 
also. While there is a possibility of variable parameter values 
resulting in ill conditioning of the system, this does not appear 
to be a real problem. 

5.0 Conclusions and Suggestions for Further Research 

The methods we have developed here appear to be very useful 
for the purposes we consider, that of approximating surfaces from 
scattered data efficiently for use in subsequent computations. 
Which of the three algorithms one might employ to obtain the 
approximation depends on several matters that are peculiar to the 
data being approximated, as well as the computational 
reguirements and resources available. If a reasonable 
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approximation is required with no heavy burden on subsequent use, 
the greedy knot selection process will probably yield a suitable 
result. If the final use imposes a high value on efficiency of 
evaluation, no doubt use of the two schemes giving optimized knot 
locations will look attractive. From our modest experiments it 
seems that use of variable parameter values at the knots adds 
approximation power beyond its cost. 

To begin with, it is desirable to carry out the 
investigation with many more sets of data. Exploration of the 
process as we have done here is very important, using known 
underlying surfaces being approximated for comparison purposes. 
However, ultimately the use of the scheme must be for 
approximation of data obtained experimentally, or from 
environmental measurements. This data is almost invariably 
subject to error. While we have does some experimentation with 
such data (e, the glacier data) , much remains to be done. 

There are a number of directions in which this research can 
be extended. One idea we have explored slightly is that of using 
some measure of smoothness of the surface, in connection with the 
rms error at the data points, to determine when to stop the knot 
addition process in the greedy scheme. A reasonable stopping 
criterion is a necessity in approximating real-world data, 
especially if the error characteristics are largely unknown. We 
have computed the approximate value of the thin plate functional 
for these surfaces with the idea that an significant increase in 
the value of the functional accompanied by only slight decrease 
in the rms error may signal that complexity is being added to the 
surface without actually improving the fit to the underlying 
surface by much. We believe it will probably be necessary to 
monitor the values over some small numbers of knots, say over 5 
or so consecutive numbers of knots. We intend to explore this as 
a potential stopping criterion. 

In certain sets of data it may desirable to include a 
smoothing term along with the rms error at the data points as 
part of the objective function in the knot location optimization 
schemes. One could take the objective function to be a convex 
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combination of the rms error and the value of some functional 
related to smoothness of the surface, such as the thin plate 
functional. There could be several reasons for this being 
desirable, but one is that if there are relative voids in the 
data, addition of a smoothing term would tend to give some 
control over the behavior of the function in such regions. There 
are many unknown factors in such a process. There are numerous 
cases where such objective functions have been found useful. 
See, for example, [HS91]. 

The particular form of the measure of smoothness probably 
depends on the application, and the use of the thin plate 
functional, while convenient and useful in many cases, may not be 
the proper one for environmental applications, for example. For 
meteorological problems it has been found that functionals 
corresponding to higher powers of the Laplacian seem to be 
appropriate [FR90]. Whatever the form of the measure of 
smoothness, the appropriate choice of weighting between the rms 
errors and the smoothness will also have to be discovered. 
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Figure Captions 



a b 

All figures are oriented (sideways) as 

c d 

Figure 1: a) The 500 data points and the subset of 50 of them 
chosen as knot points as generated by the greedy algorithm are 
shown. Minimum knot separation enforced was 0.05. Contours of 
the parent function sampled for the data are also shown. b) A 
perspective plot of the parent function, viewed from a point in 
the first quadrant. c) A perspective plot of the plot of the 
approximating multiquadric function computed as the least squares 
approximation, d) Contours of the approximating function. 

Figure 2: The surface representing two terms corresponding to 
two nearly repeated knots in a least squares approximation with 
optimized knot locations. The region is above the [0,1]^ square, 
on which the surface is sampled. b) The surface derived by 
averaging the location of the two knots and adding the 
coefficients. c) Contours of the surface corresponding to the 
two terms in part a. d) The contours of the single multiquadric 
term in part b. 

Figure 3; a) The parent surface sampled at the 500 points show 
in Figure la. b) The approximating least squares multiquadric 
with knots at the initial guess points, as in Figure la. c) The 
surface corresponding to the least squares approximation using 
the optimized knot locations. d) The initial guesses and the 
optimized knot locations. 

Figure 4; The surface representing two terms corresponding to 
two nearly repeated knots in a least squares approximation with 
optimized knot locations, each with optimized multiquadric 
parameter value; the optimized parameter values are essentially 
the same. The region is above the [0,1]^ square on which the 
surface is sampled. b) The surface derived by averaging the 
location of the two knots and adding the coefficients. c) 
Contours of the surface corresponding to the two terms in part a. 
d) The contours of the single multiquadric term in part b. 
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Figure 5: a) The 100 data points and the subset of 12 of them 
chosen as knot points as generated by the greedy algorithm are 
shown. Minimum knot separation enforced was 0.2. Contours of 
the parent function sampled for the data are also shown. b) A 
perspective plot of the parent function, viewed from a point in 
the first quadrant. c) A perspective plot of the plot of the 
approximating multiquadric function computed as the least squares 
approximation, d) Contours of the approximating function. 

Figure 6 a) The parent function which was sampled at the 100 
points shown in Figure 5a. b) The least squares approximation 
constructed from the initial guess knot points, shown in part d 
as o's. c) The multiquadric approximation constructed from the 
optimized knot locations and (single) parameter value. d) The 
initial guess (o's) and optimized (x's) knot locations. 

Figure 7: a) The parent function which was sampled at the 100 
points shown in Figure 5a. b) The least squares approximation 
constructed from the initial guess knot points, shown in part d 
as o's. c) The multiquadric approximation constructed from the 
optimized knot locations and associated parameter values. d) 
The initial guess (o's) and optimized (x's) knot locations. 
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